In [ ]:
Idiots Guide to (Open) Data Science¶
Andrew Bolster
- bolster.online
- Tweets @bolster
- GitHubs at andrewbolster
- Works at AlertLogic We're hiring Threat/Exploit Researchers (Placements also available)
- Plays at Farset Labs
- THIS NOTEBOOK IS AT present.bolster.online
- Also available in source at presentgh.bolster.online
What is Data Science?¶
- "Deriving actionable business and operational insights from multi-modal data sources"
- AKA: "Turning Numbers into Other Numbers, and occasionally pretty graphs"
Data Science in Northern Ireland¶
Vibrant Corporate Ecosystem¶
- Analytics Engines
- AquaQ Analytics
- AlertLogic😜
- Kainos
- Flexera|BDNA
- First Derivitives
- Neueda
- Sensum
- BrainwaveBank
- Many Many Many More
Open Source / Meetup Ecosystem too¶
Hint: If you really want to learn data science, go to/get involved in some of these
Hint the Second: These are the best ways to get a job in the field
Engaged local Government/ComVol support¶
NI Data Science on the world stage¶
What is Open Data?¶
Open data and content can be freely used, modified, and shared by anyone for any purpose - The Open Definition
Open data is the idea that some data should be freely available to everyone to use and republish as they wish, without restrictions from copyright, patents or other mechanisms of control. - Wikipedia
Open data is data that’s available to everyone to access, use and share. Yep, even your nan.- The Open Data Institute
Key Principles of Open Data¶
- Open
- Accessible
- Available
- Understandable
- Traceable
Open data should be easy to access. Especially for your nan.¶
Open data is only useful if it’s shared in ways that people can actually understand. It needs to be shared in a standardised format and easily traced back to where it came from.
Open data isn’t the same as big data, but big data can be open data too¶
When people talk about ‘big data’ they mean a lot of data. Obviously. But ‘big data’, like government statistics on health care, is often closed to the public (or expensive to access, like the terrifyingly-named TWITTER FIREHOSE). Opening up big data lets people use it to spot trends, fill gaps and improve services.
Oh, and it’s also not the same as ‘shared data’¶
If you’re worried about big companies being fed all your private details, that’s got nothing to do with open data. Groups sharing information with each other is different from opening it up for all to access. Your private data should only be open if you choose to share it. (But if you want to know who’s accessing or sharing your data, open data can help.)
Open data is good for democracy¶
If citizens know about their governments they can hold leaders to account, make more informed decisions and demand better services. Open data can also help governments stay on their toes and make better policies for society, the economy and the environment.
Open data can help fight crime¶
It’s helped people in London to track stolen bikes and police in Vancouver to stay one step ahead of criminals.
Open data is good for your health¶
With tools like FoodTradeMenu using it to help restaurants make sure they don’t serve you food you’re allergic to without realising.
Open data helps you get around your city, and saves you money¶
Apps like CityMapper use open data from groups like Transport for London to help you find the quickest and cheapest way to get from A to B. Even maps can be open, like OpenStreetMap, which powers map data for websites and humanitarian crises relief around the world.
And how shit is NI? Cus we're shit at everything right?¶
NOTBAD.JPG¶
What is Data Driven Journalism?¶
‘Data journalism’ only differs from ‘words journalism’ in that we use a different kit. We all sniff out, report, and relate stories for a living. It’s like ‘photo journalism’; just swap the camera for a laptop. - Brian Boyer, Chicago Tribune
[...] providing information and analysis to help inform us all about important issues of the day. - Aron Pilhofer, New York Times
Using Open Data and Data Science principles to ask, analyse and answer complex or contentious questions or areas of inquiry using available evidence. - Bolster, Here
What do you what to know?¶
This is usually the wrong question;
What do you want to know more about?¶
What relevant data is available?¶
So what do we think we have at the moment?¶
- List of schools with at least council/constitutency locations and postcodes
- Enrolements per school for 2016/17
- Demographic comparisons (R/FSM/SEN 0+ or 5)
- A 'churn' rate (newcomers/movers/immigration)
What else could we do with?¶
- More years, ideally same format
- What schools have fewer places per/k council pop?
- Correlations to deprivation
- Correlations to monoculture/political alignment
- Correlations to population change
- Correlations to achievement/mobility
- Correlations to claiment count etc.
I think we can get most of these... but first things first
Extraction, Transformation and Loading¶
- Get the data
- Clean the data
- Store the data
Gov Open Data is notoriously difficult to 'wander around'
Extraction¶
When in doubt, bring in pandas
Easy mode: read_csv¶
In [5]:
import pandas as pd
url_from_odni = "https://www.opendatani.gov.uk/dataset/6058be29-b2e1-4253-bab8-8a018568560a/resource/b4fcfed1-2dc1-4f61-8968-573535522d53/download/school-level-post-primary-reference-data.csv"
pd.read_csv(url_from_odni)
Out[5]:
Hardmode: API's¶
In [6]:
from ckanapi import RemoteCKAN # it's on pip
def build_odni_connection():
"""Be nice to OpenDataNI and tell them how old I am. (And that it's me)"""
version_no = (pd.to_datetime('now') -
pd.to_datetime('1988/05/17')).days/365
ua = f'@Bolster/{version_no:.2f} (+http://bolster.online/)'
return RemoteCKAN('https://www.opendatani.gov.uk/', user_agent=ua)
In [7]:
odni = build_odni_connection()
odni.action.package_list()
Out[7]:
In [8]:
resource = 'school-census-post-primary'
dfs = {}
for dataset in odni.action.package_show(id=resource)['resources']:
df = pd.read_csv(dataset['url'], index_col=False)
dfs[dataset['name']] = df
print(list(dfs.keys()))
In [9]:
dataset
Out[9]:
Those who ignore history...¶
There's a problem; We want to know about change over time.
Dept of Ed have this data 2009-18 on their website in Excel format
OpenDataNI do not 😭
We must deeeeel weeeth eet¶
In [10]:
test_data = "https://www.education-ni.gov.uk/sites/default/files/publications/education/School%20level%20-%20post%20primary%20schools%20data%20supp%201718.xls"
pd.read_excel(test_data).head()
Out[10]:
Common Data Extraction Hellscapes¶
- Inconsistent naming conventions
- Inconsistent format conventions
- Inconsistent schema conventions
- Inconsistent lexical conventions
- Inconsistent policy changes
- Shifting definitions/classifications
- Lack of meta-linkage during all of the above changes
Spot the difference...¶
School%20level%20-%20post%20primary%20schools%20data%20supp%201718.xlsSchool%20level%20-%20post%20primary%20schools%20data%201617%20supp_0.XLS
- order of fields (
supp 1718vs201617 supp_0) - extra junk (
suppvssupp_0) - inconstent extensions (
xlsvsXLS)
Web Crawling to the rescue!¶
In [11]:
from bs4 import BeautifulSoup
import requests
base_url="https://www.education-ni.gov.uk"
listing_path = "/articles/school-enrolments-school-level-data"
soup = BeautifulSoup(requests.get(base_url+listing_path).text, 'lxml')
for link in soup.find_all('a'):
if 'School enrolments - school level data 20' in ' '.join(map(str,link.contents)):
print(link)
Then for each of those, get the href associated with the year¶
In [12]:
soup = BeautifulSoup(requests.get(base_url+listing_path).text, 'lxml')
year_urls = {}
for link in soup.find_all('a'):
contents = ' '.join(map(str,link.contents))
if 'School enrolments - school level data 20' in contents:
year = link.get('href')[-6:-2]
year_urls[int(year)] = link.get('href')
year_urls
Out[12]:
Then walk those pages and grab the Excel files...¶
Using urllib.request.urlretrieve to store them on the local filesystem
In [13]:
import urllib.request
from pathlib import Path
from collections import defaultdict
year_files = defaultdict(list)
dest = Path('./data/education-ni/')
for year, path in year_urls.items():
dest.joinpath(str(year)).mkdir(parents=True, exist_ok=True)
soup = BeautifulSoup(requests.get(base_url+path).text, 'lxml')
for link in soup.find_all('a'):
href=link.get('href',"")
filename = href.split('/')[-1]
contents = ' '.join(map(str,link.contents))
if 'xls' in href.split('.')[-1].lower():
# See this requirement right 👆
# Mix of xls, XLSX, xlsx and XLS
if not dest.joinpath(str(year),filename).exists():
urllib.request.urlretrieve(href, dest.joinpath(str(year),filename))
year_files[year].append(filename)
print(f"{year}:{len(year_files[year])} docs")
Transformation¶
We've technically 'Extracted' the data we wanted from the source,
Clean and twist the data until it's useful
This is the most frustrating part of being a data scientist/engineer...
Quick inspection of filename 'conventions'¶
In [14]:
for year, filenames in year_files.items():
for file in filenames:
print(f"{year}:/{file}")
Zero in on just post-primary data¶
In [15]:
for year, filenames in year_files.items():
for file in filenames:
if 'post' in file:
print(f"{year}:/{file}")
Don't try and tackle it all at once¶
Right, lets limit our scope to look at places available to get a better understanding of the data
In [16]:
for year, filenames in year_files.items():
print(f"Year: {year}")
for file in filenames:
if file.lower().startswith('available'):
print(file)
Try and identify subclasses of datasets that may be tackled automagically¶
Ok, what about actually parsing the excel sheets?
Will they be consistent?
In [17]:
sheets = defaultdict(list)
for year, files in year_files.items():
for file in files:
if 'post' in file:
xls = pd.ExcelFile(f'data/education-ni/{year}/{file}')
sheets[year].extend(xls.sheet_names)
dict(sheets)
Out[17]:
In [18]:
from collections import Counter
all_sheet_names = Counter([_ for d in sheets.values() for _ in d])
all_sheet_names.most_common()
Out[18]:
In [19]:
df = pd.DataFrame.from_dict({
year: [sn in sheet_names for sn in all_sheet_names]
for year, sheet_names in sheets.items()
}, orient='index')
df.columns=all_sheet_names.keys()
df.T
Out[19]:
Visualisation isn't just for output data...¶
Use visualisation tools to validate your transformation assumptions...
All. The. Time.
In [20]:
import seaborn as sns
f,ax = plt.subplots(figsize=graph_figsize)
sns.heatmap(df.T, ax=ax)
Out[20]:
Sometimes understanding the data is more painful (and ugly) than expected¶
This was the only document that seemed to explain what 'supernumerary' actually means
It means 'extra'...
Using multiple tools simultanously can be problematic...¶
Excel on Mac helpfully puts down a lock-entry on open excel files which prevent it from being opened by other things.
Pandas occasionally has unhelpful error messages in this regard
Sometimes you have to give up on stuff because it's too dangerous an area to walk in to.¶
Special Educational Needs classification policy is >>100pg policy that appears to change every year.
I am not an expert on education policy so I'm not sticking my foot in that one without being sure I understand the implications.
It's your responsibility as a Data Science to be able to make that call, and know your limits.
"Machine Readable"¶
<img src="img/hashs.png", width="90%"/>
In [21]:
def parse_reference_table(xls):
"""
From an ExcelFile, clean up:
* School Mgmt Type disaster
* Inconsistent header depth
* Multi-row header names
* inconsistent headers (ref_key_map)
* inconsistent col order
* inconsistent caps/spacing (strip|lower)
"""
cols= [
'de ref',
'school name',
'school type',
'address 1',
'postcode',
'urban_rural',
'school management type',
'district council',
'parliamentary constituency',
'town'
]
categories = [
'school type',
'urban_rural',
'school management type',
'district council',
'parliamentary constituency',
'town'
]
ref_key_map={
'denino':'de ref',
'urban/ rural': 'urban_rural',
'schoolname': 'school name'
}
reference_value_rename = {
'school management type':{
'gmi':'integrated',
'controlled integrated':'integrated',
'roman catholic maintained':'rc maintained',
'grant maintained integrated':'integrated',
'voluntary - other managed':'voluntary',
'voluntary - rc managed':'voluntary',
'catholic maintained':'rc maintained'
}
}
join_n_strip_n_lower = lambda l: ' '.join(l).strip().lower()
if 'reference data' in xls.sheet_names:
df = pd.read_excel(xls, 'reference data', header=None)
h_range = 2 if isinstance(df.ix[3,0], int) else 3
try:
df.columns=df.ix[1:h_range].fillna('').apply(join_n_strip_n_lower, axis=0).values
df.rename(columns=ref_key_map, inplace=True)
df = df.drop(df.index[0:h_range+1]).reset_index(drop=True)
df = df[cols]
df['de ref'] = df['de ref'].astype(int)
df.set_index('de ref', inplace=True)
for c in df:
df[c]=df[c].str.lower().str.strip()
df.replace(reference_value_rename, inplace=True)
for c in categories:
df[c].fillna('NA', inplace=True)
df[c] = df[c].astype('category')
except TypeError as e:
print(e)
else:
df=None
return df
def parse_enrolments_table(xls):
"""From an ExcelFile, clean up:
* Inconsistent header depth
* fucked up nans/nulls all over the place
* inconsistent *footer* depth...
* Set de ref as join index and drop pointless fields
"""
def join_n_strip_n_lower(l): return ' '.join(l).strip().lower()
def strip_n_lower(s): return s.strip().lower()
def unyearify(s): return int(s.replace('year ', ''))
if 'enrolments' in xls.sheet_names:
df = pd.read_excel(xls, 'enrolments', header=None, skip_footer=5)
h_range = 2 if isinstance(df.ix[3, 0], int) else 3
try:
df.columns = df.ix[3].fillna('').apply(strip_n_lower).values
df = df.drop(df.index[0:h_range+1]).reset_index(drop=True)
df.dropna(how='all', inplace=True, axis=0)
df['de ref'] = df['de ref'].astype(int)
df.drop('schoolname', axis=1, inplace=True)
df.drop('total pupils', axis=1, inplace=True)
df.set_index('de ref', inplace=True)
df.rename(columns=unyearify, inplace=True)
df = df.astype(float)
except TypeError as e:
print(e)
else:
df = None
return df
def parse_fsm_table(xls):
"""From an ExcelFile, clean up:
* Inconsistent header depth
* fucked up nans/nulls all over the place
* inconsistent *footer* depth...
* Set de ref as join index and drop pointless fields
"""
join_n_strip_n_lower = lambda l: ' '.join(l).strip().lower()
strip_n_lower = lambda s: s.strip().lower()
unyearify = lambda s: int(s.replace('year ',''))
if 'free school meals' in xls.sheet_names:
df = pd.read_excel(xls, 'free school meals', header=None, skip_footer=5)
h_range = 2 if isinstance(df.ix[3,0], int) else 3
try:
df.columns=df.ix[3].fillna('').apply(strip_n_lower).values
df = df.drop(df.index[0:h_range+1]).reset_index(drop=True)
df.dropna(how='all', inplace=True, axis=0)
df['de ref'] = df['de ref'].astype(int)
df.drop('schoolname',axis=1, inplace=True)
df.drop('free school meals', axis=1, inplace=True)
df.set_index('de ref', inplace=True)
df.replace('#',pd.np.nan, inplace=True) # # = Undisclosed
df.replace('*',2.0, inplace=True) # * == < 5
df.replace('!',1, inplace=True) # ! avoid identification, so it's prob one or two
df=df.astype(float)
except TypeError as e:
print(e)
else:
df=None
return df
def parse_available_table(xls):
"""From an ExcelFile, clean up:
* Inconsistent header depth
* fucked up nans/nulls all over the place
* inconsistent *footer* depth...
* Set de ref as join index and drop pointless fields
* Totally different schemas between years
* Inconsistent metric naming
* non numerical data flags (*/!)
"""
ref_key_map={
'schoolname': 'school name',
'total unfilled places': 'available places',
'unfilled places': 'available places',
'total approved enrolment number': 'approved enrolments'
}
join_n_strip_n_lower = lambda l: ' '.join(l).strip().lower()
if 'School level data' in xls.sheet_names:
df = pd.read_excel(xls, 'School level data', header=None)
h_range = 2 if isinstance(df.ix[3,0], int) else 3
elif 'unfilled places' in xls.sheet_names:
df = pd.read_excel(xls, 'unfilled places', header=None)
h_range = 2 if isinstance(df.ix[3,0], int) else 3
else:
df=None
if df is not None:
try:
df.columns=df.ix[1:h_range].fillna('').apply(join_n_strip_n_lower, axis=0).values
df.rename(columns=ref_key_map, inplace=True)
df = df.drop(df.index[0:h_range+1]).reset_index(drop=True)
df=df.applymap(lambda x: np.nan if isinstance(x, str) and x.isspace() else x)
df.dropna(how='all', axis=1, inplace=True)
df.dropna(how='any', axis=0, inplace=True)
if df.shape[1] == 6: # recent doesn't have fecking headers
cols = list(df.columns)
cols[0] = 'de ref'
cols[1] = 'school name'
df.columns=cols
df.drop('school name', axis=1, inplace=True)
df['de ref'] = df['de ref'].astype(int)
df.set_index('de ref', inplace=True)
df.replace('*',2.0, inplace=True) # * == < 5
df.replace('!',1, inplace=True) # ! avoid identification, so it's prob one or two
df.dropna(how='all', inplace=True, axis=1)
df.astype(int, inplace=True)
except TypeError as e:
print(e)
return df
In [22]:
re_dfs={}
av_dfs={}
en_dfs={}
fsm_dfs={}
for year, files in year_files.items():
for file in files:
if 'post' in file:
xls = pd.ExcelFile(f'data/education-ni/{year}/{file}')
df = parse_reference_table(xls)
if df is not None:
print(f'Got reference data for {year}')
re_dfs[year]=df
df = parse_enrolments_table(xls)
if df is not None:
print(f'Got enrolment data for {year}')
en_dfs[year]=df
df = parse_available_table(xls)
if df is not None:
print(f'Got available data for {year}')
av_dfs[year]=df
df = parse_fsm_table(xls)
if df is not None:
print(f'Got fsm data for {year}')
fsm_dfs[year]=df
reference = pd.Panel(re_dfs).sort_index()
available = pd.Panel(av_dfs).sort_index()
enrolment = pd.Panel(en_dfs).sort_index()
fsm = pd.Panel(fsm_dfs).sort_index()
Bring it all together¶
Panels are 3D Dataframes, i.e. [Year x School x Metric]
- Reference: School name, address, denomination, type, etc.
- Available: Allocated places information
- Enrolment: Actual student numbers
- FSM: Free School Meals: Laziest possible deprivation measure, probably doesn't need to be a panel
In [23]:
enrolment.sum()
Out[23]:
In [24]:
enrolment.sum().sum().iplot()
In [25]:
enrolment.sum().iplot(kind='bar')
In [26]:
enrolment.sum().T.iplot(kind='bar',barmode='stack')
In [28]:
# What is wrong with this "sum" table?
enrolment.sum(axis=2).head()
Out[28]:
In [29]:
[10]+[np.nan]*10
Out[29]:
In [30]:
sum([10]+[np.nan]*10)
Out[30]:
In [31]:
pd.Series([10]+[np.nan]*10).sum()
Out[31]:
In [32]:
pd.Series([10]+[np.nan]*10).mean()
Out[32]:
In [33]:
pd.Series([np.nan]*10).sum()
Out[33]:
In [34]:
pd.Series([10]+[0.0]*10).mean()
Out[34]:
In [35]:
# These 0's are actually sum([np.nan])
# The school was closed in 2014
enrolment.sum(axis=2).head()
Out[35]:
In [36]:
reference[2013].loc[1210266]
Out[36]:
In [37]:
reference[2014].loc[1210266]
Out[37]:
In [38]:
# It *should* look like this
enrolment.sum(axis=2).head().replace(0.0, np.nan)
Out[38]:
AUDIENCE PARTICIPATION¶
Which of the below is the correct answer to: "What is the average number of pupils per school in Northern Ireland from 2009 to 2017?"
In [39]:
#A Sum of years, swap 0 for nan, mean over schools, mean over classes
enrolment.sum(axis=1).replace(0.0, np.nan).mean(axis=1).mean()
Out[39]:
In [40]:
#B Sum over classes, mean over schools, swap 0 for nan, mean over years
enrolment.sum(axis=2).mean().replace(0.0,np.nan).mean()
Out[40]:
In [41]:
#C Sum over classes, swap 0.0 for nan, mean over schools, mean over years
enrolment.sum(axis=2).replace(0.0, np.nan).mean().mean()
Out[41]:
C: Sum over classes, swap 0.0 for nan, mean over schools, mean over years!¶
But why not B, it's doing the same thing right?
In [42]:
enrolment.sum(axis=2).head()
Out[42]:
In [43]:
enrolment.sum(axis=2).head().mean().to_frame().T
Out[43]:
In [44]:
enrolment.sum(axis=2).replace(0.0,np.nan).head()
Out[44]:
In [45]:
enrolment.sum(axis=2).replace(0.0,np.nan).mean().to_frame().T
Out[45]:
It's your responsibility as a Data Scientist to understand what questions your algorithms are really answering, and make sure that the people who may be answering those question are crystal clear about what your data means.
Data Science is as much about comprehension and communication as it is about statistics and programming.
Return of regular schedule - What happened to average cohort sizes in 2015?¶
In [46]:
enrolment.mean(axis=2).mean().iplot(title='Annual Average Cohort (school/class/year)')
In [47]:
enrolment.mean(axis=2).T.iplot(title='Annual Average Cohort (school/class/year) per school', legend=False)
In [48]:
annual_school_enrolment = enrolment.sum(axis=2).T
annual_school_enrolment.head()
Out[48]:
In [109]:
name_lookup = reference.minor_xs('school name').T.describe().ix['top']
annual_school_enrolment.rename(columns=name_lookup).iplot(legend=False)
In [50]:
first = annual_school_enrolment.replace(0.0,np.nan).apply(pd.Series.first_valid_index).replace(2009,np.nan)
last = annual_school_enrolment.replace(0.0,np.nan).apply(pd.Series.last_valid_index).replace(2017,np.nan)
first.dropna()
Out[50]:
In [51]:
close_vec=pd.DataFrame({
'opened':first.groupby(first).size(),
'closed':last.groupby(last).size()}
)
close_vec
Out[51]:
In [52]:
close_vec.cumsum().fillna(method='ffill').iplot(title="Cumulative number of schools opened/closed since 2009")
Skippable Aside: Slicing Panels¶
In [53]:
up_to_ks3=enrolment.minor_axis<=10
up_to_ks3
Out[53]:
In [54]:
enrolment[:,:,up_to_ks3].sum()
Out[54]:
In [55]:
cowtown=reference.major_axis[(reference.minor_xs('town')=='ballymoney').any(axis=1)]
cowtown
Out[55]:
In [56]:
enrolment.sum(axis=2).loc[cowtown].sum()
Out[56]:
Mapping¶
- Maps are hard in NI
- No Open Data mapping between Postcode -> Boundary
- Also doesn't help that ESRI's API's have been broken recently....
However... we do get towns and administrative boundaries... from something called the "Gazetter"
In [57]:
if not Path('data/gztr').exists():
gztr_json=requests.get("http://osni-spatial-ni.opendata.arcgis.com/datasets/117e5c3d0f0b41208d1caec1ddbd1330_1.geojson").json()
gztr = gp.GeoDataFrame.from_features(gztr_json).set_index('OBJECTID')
gztr['PLACENAME'] = gztr['PLACENAME'].str.lower()
gztr.to_file('data/gztr')
gztr = gp.read_file('data/gztr')
gztr.head()
Out[57]:
In [58]:
gztr.plot()
Out[58]:
And we have a base map...¶
In [59]:
base_map_zip="http://osni-spatial-ni.opendata.arcgis.com/datasets/d9dfdaf77847401e81efc9471dcd09e1_0.zip"
base_map_shp = "OSNI_Open_Data_Largescale_Boundaries__NI_Outline.shp"
if not Path('data/'+base_map_shp).exists():
urllib.request.urlretrieve(base_map_zip, 'data/_tmp.zip')
with zipfile.ZipFile('data/_tmp.zip') as z:
z.extractall('data/')
Path('data/_tmp.zip').unlink()
base_map=gp.GeoDataFrame.from_file('data/'+base_map_shp)
In [60]:
f,ax = plt.subplots(figsize=graph_figsize)
base_map.plot(ax=ax, color='gray')
gztr.plot(ax=ax, color='b')
Out[60]:
In [61]:
from scipy.stats import mode
school_towns = reference.minor_xs('town').mode(axis=1)[0].value_counts()
school_towns
Out[61]:
In [62]:
school_towns_gp = gztr[gztr['PLACENAME'].isin(school_towns.keys())].set_index('PLACENAME')
school_towns_gp.head()
Out[62]:
In [63]:
school_towns_gp['school_count'] = school_towns.astype(float)
school_towns_gp.head()
Out[63]:
In [64]:
f,ax = plt.subplots(figsize=graph_figsize)
base_map.plot(ax=ax, color='gray')
school_towns_gp.plot(ax=ax, column='school_count', legend=True)
Out[64]:
Points aren't really very good for seeing trends though...¶
Need some regions
In [65]:
cons_map_zip="http://osni-spatial-ni.opendata.arcgis.com/datasets/563dc2ec3d9943428e3fe68966d40deb_3.zip"
cons_map_shp = "OSNI_Open_Data_Largescale_Boundaries__Parliamentary_Constituencies_2008.shp"
if not Path('data/'+cons_map_shp).exists():
urllib.request.urlretrieve(cons_map_zip, 'data/_tmp.zip')
with zipfile.ZipFile('data/_tmp.zip') as z:
z.extractall('data/')
Path('data/_tmp.zip').unlink()
cons_map=gp.GeoDataFrame.from_file('data/'+cons_map_shp)
cons_map.rename(columns={'PC_NAME':'constituency'}, inplace=True)
cons_map.drop(['OBJECTID','PC_ID'], axis=1, inplace=True)
cons_map['constituency'] = cons_map['constituency'].str.lower().str.strip()
cons_map.set_index('constituency', inplace=True)
cons_map.head()
Out[65]:
In [66]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_map.plot(ax=ax)
cons_map.apply(lambda x: ax.annotate(x.name,
xy=x.geometry.centroid.coords[0],
ha='center'),
axis=1);
In [67]:
reference.minor_xs('parliamentary constituency').head()
Out[67]:
In [68]:
cons = set(reference.minor_xs('parliamentary constituency').values.ravel('k'))
en_df= pd.DataFrame.from_dict({
con:
enrolment[:,reference.major_axis[(reference.minor_xs('parliamentary constituency')==con).any(axis=1)],:].sum().sum()
for con in cons
}).T.drop(np.nan)
av_df= pd.DataFrame.from_dict({
con:
available[:,reference.major_axis[(reference.minor_xs('parliamentary constituency')==con).any(axis=1)],'available places'].sum()
for con in cons
}).T.drop(np.nan)
cons_stats = gp.GeoDataFrame(geometry=cons_map.geometry)
for c in av_df:
cons_stats[f"av_{c}"] = av_df[c]
cons_stats[f"av_{c}_rat"] = (av_df[c]/en_df[c])
for c in en_df:
cons_stats[f"en_{c}"] = en_df[c]
cons_stats[f"en_{c}_pk"] = en_df[c]/cons_map["Area_sqkm"]
In [69]:
cons_stats.head()
Out[69]:
In [70]:
f,ax = plt.subplots(figsize=graph_figsize)
ax.set_title("Per Constituency Total Enrolments [2017]")
cons_stats.plot(column='en_2017',ax=ax, legend=True)
cons_stats.apply(lambda x:
ax.annotate(x.name,
xy=x.geometry.centroid.coords[0],
ha='center',
color='#eeeeee',
fontsize=50*x.geometry.area
),
axis=1);
In [71]:
def annotate_cons(x, ax=ax):
ax.annotate(x.name,
xy=x.geometry.centroid.coords[0],
ha='center',
color='#eeeeee',
fontsize=50*x.geometry.area
)
In [72]:
f,ax = plt.subplots(figsize=graph_figsize)
ax.set_title("Per Constituency Total Enrolments per $km^2$[2017]")
cons_stats.plot(column='en_2017_pk',ax=ax, legend=True)
cons_stats.apply(annotate_cons, ax=ax, axis=1);
In [73]:
cons_stats["en_2017_pk"].sort_values(ascending=False)
Out[73]:
But what does it all mean¶
Gather data from other sources to put things in context;
How many people are in each region?
In [74]:
if not Path('data/cons_pop.csv').exists():
for dataset in odni.action.package_show(id='population-estimates-for-northern-ireland')['resources']:
if dataset['name'] == "Parliamentary Constituencies by single year of age and gender (mid-2001 to mid-2017)":
cons_pop = pd.read_csv(dataset['url'])
cons_pop.to_csv("data/cons_pop.csv", index=False)
cons_pop.head()
cons_pop = pd.read_csv('data/cons_pop.csv')
cons_pop['Mid_Year_Ending'] = cons_pop.Mid_Year_Ending.astype(int)
cons_pop['Population_Estimate'] = cons_pop.Population_Estimate.astype(float)
cons_pop['Age'] = cons_pop.Age.astype(int)
cons_pop.rename(columns={'Geo_Name':'constituency'}, inplace=True)
cons_pop['constituency']= cons_pop.constituency.str.strip().str.lower()
In [75]:
cons_pop[(cons_pop.Gender == 'All Persons') & (cons_pop.Mid_Year_Ending == 2016)].head()
Out[75]:
What's the distribution of ages in NI?¶
In [76]:
f,ax = plt.subplots(figsize=graph_figsize)
ax.set_title('Population Distribution (2016)')
cons_pop[(cons_pop.Gender == 'All Persons') & (cons_pop.Mid_Year_Ending == 2016)]\
.groupby('Age')['Population_Estimate'].sum().plot.line(ax=ax)
Out[76]:
Constituency breakdown¶
In [77]:
f,ax = plt.subplots(figsize=graph_figsize)
ax.set_title('Constituency Population Distribution (2016)')
cons_pop[(cons_pop.Gender == 'All Persons') & (cons_pop.Mid_Year_Ending == 2016)]\
.groupby(['Age','constituency'])['Population_Estimate'].sum().unstack().plot.bar(stacked=True, ax=ax)
Out[77]:
That's not a great way of summarising age distributions....¶
In [78]:
def idxquantile(s, q=0.5, *args, **kwargs):
qv = s.quantile(q, *args, **kwargs)
return (s.sort_values()[::-1] <= qv).idxmax()
f,ax = plt.subplots(figsize=graph_figsize)
ax.set_title('Cumulative Population Distribution (2016)')
pop_cumsum = cons_pop[(cons_pop.Gender == 'All Persons') & (cons_pop.Mid_Year_Ending == 2016)]\
.groupby('Age')['Population_Estimate'].sum().cumsum()
pop_cumsum.plot(ax=ax)
age_quartiles = [idxquantile(pop_cumsum, q) for q in [0.25,0.5,0.75]]
ax.vlines(age_quartiles, ymin=pop_cumsum.min(), ymax = pop_cumsum[age_quartiles])
age_quartiles
Out[78]:
In [79]:
f,ax = plt.subplots(figsize=graph_figsize)
ax.set_title('Cumulative Population Distribution (2009-2016)')
yrs=sorted(cons_pop.Mid_Year_Ending.unique())
colors = plt.cm.jet(np.linspace(0,1,len(yrs)))
for yr,c in zip(yrs,colors):
pop_cumsum = cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.Mid_Year_Ending == yr)]\
.groupby('Age')['Population_Estimate'].sum().cumsum()
pop_cumsum.plot(ax=ax, label=yr, c=c)
age_quartiles = [idxquantile(pop_cumsum, q) for q in [0.25,0.5,0.75]]
ax.vlines(age_quartiles, ymin=pop_cumsum.min(), ymax = pop_cumsum[age_quartiles],color=c,alpha=0.2)
ax.legend()
Out[79]:
In [80]:
def constituency_age_profile(constituency=None):
# I have a hunch we might was to do this
# with lists....
if constituency is None:
constituency = cons_pop.constituency.unique()
label = f"Population profile for NI"
elif not isinstance(constituency, list):
label = f"Population profile for {constituency}"
constituency = [constituency]
else:
label = f"Population profile for {constituency}"
f,ax = plt.subplots(figsize=graph_figsize)
yrs=sorted(cons_pop.Mid_Year_Ending.unique())
colors = plt.cm.jet(np.linspace(0,1,len(yrs)))
for yr,c in zip(yrs,colors):
pop_cumsum = cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.constituency.isin(constituency))\
& (cons_pop.Mid_Year_Ending == yr)]\
.groupby('Age')['Population_Estimate'].sum().cumsum()
pop_cumsum.plot(ax=ax, label=yr, c=c)
age_quartiles = [idxquantile(pop_cumsum, q) for q in [0.25,0.5,0.75]]
ax.vlines(age_quartiles, ymin=pop_cumsum.min(), ymax = pop_cumsum[age_quartiles],color=c,alpha=0.2)
for i,age in enumerate(age_quartiles):
ax.annotate(f"<={age}", xy=(age, i))
ax.legend()
ax.set_title(label)
South Belfast: Kinky Constituency¶
In [81]:
constituency_age_profile('belfast south')
Newry and Armagh: Middle Age Mourneing¶
In [82]:
constituency_age_profile('newry and armagh')
Foyled Again: Cross-generational Brain Drain¶
In [83]:
constituency_age_profile('foyle')
North Antrim: 40 is the new 30¶
In [84]:
constituency_age_profile('north antrim')
Growing Up in the (Lagan) Valley¶
In [85]:
constituency_age_profile('lagan valley')
Fast times at Belfast Unis/Colleges¶
In [86]:
df = cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.Mid_Year_Ending == 2016)]\
.groupby(['Age','constituency'])['Population_Estimate'].sum().unstack()
df.iplot()
In [87]:
f,ax = plt.subplots(figsize=graph_figsize)
sns.heatmap(df.T, ax=ax)
Out[87]:
In [88]:
def colnorm(df):
return (df-df.min())/(df.max()-df.min())
colnorm(df).iplot()
In [89]:
f,ax = plt.subplots(figsize=graph_figsize)
sns.heatmap(colnorm(df).T, ax=ax)
Out[89]:
In [90]:
def year_age_profile(yr=2016):
constituencies = list(cons_pop.groupby("constituency")['Population_Estimate'].sum().sort_values(ascending=True).keys())
label = f"Population profile for NI in {yr}"
f,ax = plt.subplots(figsize=graph_figsize)
colors = plt.cm.jet(np.linspace(0,1,len(constituencies)))
for constituency,c in zip(constituencies,colors):
pop_cumsum = cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.constituency==constituency)\
& (cons_pop.Mid_Year_Ending == yr)]\
.groupby('Age')['Population_Estimate'].sum().cumsum()
pop_cumsum.plot(ax=ax, label=constituency, c=c)
age_quartiles = [idxquantile(pop_cumsum, q) for q in [0.25,0.5,0.75]]
ax.vlines(age_quartiles, ymin=pop_cumsum.min(), ymax = pop_cumsum[age_quartiles],color=c,alpha=0.2)
ax.legend()
ax.set_title(label)
In [91]:
year_age_profile()
In [92]:
cons_pop_qilted=pd.DataFrame.from_dict(
{
qtile:cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.Mid_Year_Ending == 2016) \
& (lotile<cons_pop.Age) &(cons_pop.Age<=qtile)
].groupby('constituency')['Population_Estimate'].sum()
for lotile, qtile in zip([-1]+age_quartiles, age_quartiles+[90])
}
)
In [93]:
cons_pop_qilted.idxmax(axis=1).sort_values()
Out[93]:
In [94]:
pop_years = cons_pop.Mid_Year_Ending.unique()
for c in en_df:
if c not in pop_years:
yr = pop_years[np.abs(pop_years-c).argmin()]
else:
yr=c
_cons_pop = cons_pop[(cons_pop.Gender != 'All Persons') \
& (cons_pop.Mid_Year_Ending == yr) \
].groupby(['constituency', 'Gender'])['Population_Estimate'].sum().unstack().rename(columns=lambda c: c.lower())
_cons_pop['total'] = _cons_pop.sum(axis=1)
_cons_pop['m_per_f'] = _cons_pop['males']/_cons_pop['females']
cons_pop_qilted=pd.DataFrame.from_dict(
{
qtile:cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.Mid_Year_Ending == yr) \
& (lotile<cons_pop.Age) &(cons_pop.Age<=qtile)
].groupby('constituency')['Population_Estimate'].sum()
for lotile, qtile in zip([-1]+age_quartiles, age_quartiles+[90])
}
)
cons_stats[f"pop_{c}"] = _cons_pop['total']
cons_stats[f"pop_{c}_males"] = _cons_pop['males']
cons_stats[f"pop_{c}_females"] = _cons_pop['females']
cons_stats[f"pop_{c}_m_per_f"] = _cons_pop['m_per_f']
cons_stats[f"topqt_{c}"]=cons_pop_qilted.idxmax(axis=1)
cons_stats[f"en_{c}_pc"]=cons_stats[f"en_{c}"]/cons_pop_qilted.sum(axis=1)
cons_pop_mean=cons_pop[(cons_pop.Gender.isin(["All Persons"])) & (cons_pop.Mid_Year_Ending == yr)]
cons_pop_mean['popprod'] = cons_pop_mean[['Age','Population_Estimate']].product(axis=1)
cons_pop_mean=cons_pop_mean.groupby('constituency')[['popprod','Population_Estimate']].sum(axis=0)
cons_stats[f'age_{c}_avg'] = cons_pop_mean['popprod']/cons_pop_mean['Population_Estimate']
try:
cons_stats[f"av_{c}_pc"]=cons_stats[f"av_{c}"]/cons_pop_qilted.sum(axis=1)
except:
print(f"No data for av in {c}")
Average Age¶
In [95]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_stats.plot(column='age_2016_avg', legend=True, ax=ax)
Out[95]:
Top Quartile¶
Who rules the roost
In [96]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_stats.plot(column='topqt_2016', legend=True, ax=ax)
Out[96]:
Mars and Venus¶
In [97]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_stats.plot(column='pop_2016_m_per_f', legend=True, ax=ax)
Out[97]:
Enroled Pupils per Capita¶
Lagan Valley and South Antrim popping up again!
In [98]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_stats.plot(column='en_2016_pc', legend=True, ax=ax)
Out[98]:
Available Places per Capita¶
In [99]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_stats.plot(column='av_2016_pc', legend=True, ax=ax)
Out[99]:
Free School Meals¶
In [100]:
fsm_df= pd.DataFrame.from_dict({
con:
fsm[:,reference.major_axis[(reference.minor_xs('parliamentary constituency')==con).any(axis=1)],'fsme'].sum()
for con in cons
}).T.drop(np.nan)
for c in fsm_df:
cons_stats[f"fsm_{c}"] = fsm_df[c]
cons_stats[f"fsm_{c}_rat"] = (fsm_df[c]/en_df[c])
FSM / Enrolled¶
In [101]:
f,ax = plt.subplots(figsize=graph_figsize)
cons_stats.plot(column='fsm_2016_rat', legend=True, ax=ax)
Out[101]:
FSM is (loosly) a young persons problem¶
In [102]:
unyearify = lambda c:'_'.join([c.split('_')[0]]+c.split('_')[2:])
cons_stats_2016 = cons_stats[[c for c in cons_stats.columns if '2016' in c]].rename(columns=unyearify)
sns.pairplot(data=cons_stats_2016, hue='topqt',
vars=['age_avg','fsm_rat','pop_m_per_f'])
Out[102]:
"I'm bored, make a GIF!"¶
In [133]:
cons_stats['en_2017_pk'].sort_values()
Out[133]:
In [103]:
import matplotlib.animation as animation
unyearify = lambda c:'_'.join([c.split('_')[0]]+c.split('_')[2:])
years = list(range(2009,2017))
explanations={
'age_avg': "Average Age",
'av_pc': "Available School Places per capita",
'av_rat': "Ratio of Available Places to Actualised Enrolment",
'av': "Available School Places",
'en_pc': "Enrolled Pupils per capita",
'en_pk': "Enrolled Pupils per $km^2$",
'en': "Enrolled Pupils",
'fsm_rat': "Ratio of FSM pupils to total pupils",
'fsm': "FSM pupils",
'pop': "Population",
'pop_females': "Population (F)",
'pop_males': "Population (M)",
'pop_m_per_f': "Ratio of Males to Females",
'topqt': "Most populous age quartile (Age or younger)"
}
In [104]:
for metric in set(map(unyearify, cons_stats.columns)):
for year in years:
try:
f, ax = plt.subplots(figsize=graph_figsize)
if metric in explanations:
ax.set_title(f"{explanations[metric]}:({year})")
else:
ax.set_title(f"{metric}:({year})")
cons_stats[[c for c in cons_stats.columns if str(year) in c]+['geometry']]\
.rename(columns=unyearify).plot(column=metric, legend=True, ax=ax)
ax.axis('off')
f.tight_layout()
p = Path(f"outputs/{metric}/{year}.png")
p.parent.mkdir(parents=True, exist_ok=True)
f.savefig(p)
plt.close(f)
fps = 12
except:
print(f"{metric}:{year} failed")
finally:
if f:
plt.close(f)
In [105]:
import moviepy.editor as mpy
import moviepy.video as mpyv
clips = {}
for p in Path('outputs/').iterdir():
if not p.is_dir():
continue
metric = p.parts[-1]
file_list = [str(_p) for _p in sorted(filter(lambda s:s.suffix=='.png', p.iterdir()))]
clip = mpy.ImageSequenceClip(file_list, fps=1.5)
clip = mpyv.fx.all.freeze(clip, t='end', freeze_duration=2)
clip.write_gif(f"outputs/{metric}.gif")
clips[metric]=clip
We're getting older, or at least living longer¶
But younger people are taking over the west¶
💕 Go West for Men, or West Belfast for Women 💕¶
10% fewer men than women in West Belfast, ~2% more men than women out west
Free School Meals are rising with population¶
But the divide is getting wider per capita¶
There aren't necessarily fewer available school places either...¶
There just aren't any where people want to send their kids...
Could possibly indicate over-provisioning in yellow areas...
Population distribution isn't really changing at all¶
And as such, noones sending their kids anywhere differently despite closures...¶
🎶All together now🎶¶
In [106]:
clips.keys()
Out[106]:
In [107]:
def batch(iterable, n=1):
"""
Split an iterable into n-length batches
:param iterable:
:param n:
:return:
"""
parent_length = len(iterable)
for ndx in range(0, parent_length, n):
yield iterable[ndx:min(ndx + n, parent_length)]
ex_list = ['pop_males','pop_females']
all_gifs = [v for k,v in clips.items() if k not in ex_list]
mpy.clips_array(list(batch(all_gifs,4))).write_gif('outputs/all.gif')
Population Distribution Shifts¶
In [113]:
metric='Population Distribution'
for year in years:
try:
f, ax = plt.subplots(figsize=graph_figsize)
if metric in explanations:
ax.set_title(f"{explanations[metric]}:({year})")
else:
ax.set_title(f"{metric}:({year})")
df = cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.Mid_Year_Ending == year)]\
.groupby(['Age','constituency'])['Population_Estimate']\
.sum().unstack()
sns.heatmap(df.T, ax=ax)
f.tight_layout()
p = Path(f"outputs/{metric}/{year}.png")
p.parent.mkdir(parents=True, exist_ok=True)
f.savefig(p)
plt.close(f)
fps = 12
except:
print(f"{metric}:{year} failed")
finally:
if f:
plt.close(f)
metric='Population Distribution _normed_'
for year in years:
try:
f, ax = plt.subplots(figsize=graph_figsize)
if metric in explanations:
ax.set_title(f"{explanations[metric]}:({year})")
else:
ax.set_title(f"{metric}:({year})")
df = cons_pop[(cons_pop.Gender == 'All Persons') \
& (cons_pop.Mid_Year_Ending == year)]\
.groupby(['Age','constituency'])['Population_Estimate']\
.sum().unstack()
sns.heatmap(colnorm(df.T), ax=ax)
f.tight_layout()
p = Path(f"outputs/{metric}/{year}.png")
p.parent.mkdir(parents=True, exist_ok=True)
f.savefig(p)
plt.close(f)
fps = 12
except:
print(f"{metric}:{year} failed")
finally:
if f:
plt.close(f)
clips = {}
for p in Path('outputs/').iterdir():
if not p.is_dir():
continue
metric = p.parts[-1]
file_list = [str(_p) for _p in sorted(filter(lambda s:s.suffix=='.png', p.iterdir()))]
clip = mpy.ImageSequenceClip(file_list, fps=1.5)
clip = mpyv.fx.all.freeze(clip, t='end', freeze_duration=2)
clip.write_gif(f"outputs/{metric}.gif")
clips[metric]=clip
Conclusions¶
Hypotheses / conjectures¶
- East Antrim Constituency is literally 'little britain'
- None wants to go to school in West Belfast for some reason
- Everyone wants to go to school in North Down
- The 'places' scandal is either over blown, or more related to 'preference'
- There is a demographic wave of <11's that is gonna hit education system in 2 years. Good luck.
- South Belfast steals all-the-18-year-olds, keeps em for a decade or so and then the go home
- Foyle has been experiencing a sustained brain drain of < 35's that is now impacting services
Final Comments¶
- Data is hard
- Open Data is hard without a budget
- 80/20 rule is more like 90/10
Idiots Guide to (Open) Data Science - Thank you¶
Andrew Bolster
- bolster.online
- Tweets @bolster
- GitHubs at andrewbolster
- Works at AlertLogic We're hiring Threat/Exploit Researchers (Placements also available)
- Plays at Farset Labs
- THIS NOTEBOOK IS AT present.bolster.online
- Also available in source at presentgh.bolster.online
NOTES FOR LATER¶
Fucked Notes:¶
DoE¶
- inconsistent file names
- inconsistent extensions
- inconsistent cases in extentions 🤦
- changing schema almost every year
- recent years don't have headers for de ref/school name
- multi row headers
- Changing "All People" to "All Persons" silently...
Housing Stock Data¶
https://www.opendatani.gov.uk/dataset/ni-housing-stock
SEN scaling is complex AF¶
https://www.education-ni.gov.uk/sites/default/files/publications/de/the-code-of-practice.pdf
~ hundred page doc to work out that 5 is 'confirmed and fully stated special educational needs. Everthing less is 'suspected and in need of advice and support'
"Controlled and Maintained schools have [a budget] for this" "Grant maintained integrated schools also receive this through their annual budget"
In [ ]: